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ABSTRACT 

Resent observations and theoretical interpretations suggest that IMBHs (intermediate-mass black 
hole) are formed in the centers of young and compact star clusters born close to the center of their 
parent galaxy. Such a star cluster would sink toward the center of the galaxy, and at the same time 
stars are stripped out of the cluster by the tidal field of the parent galaxy. We investigated the orbital 
evolution of the IMBH, after its parent cluster is completely disrupted by the tidal field of the parent 
galaxy, by means of large-scale TV-body simulations. We constructed a model of the central region of 
our galaxy, with an SMBH (supermassive black hole) and Bahcall-Wolf stellar cusp, and placed an 
IMBH in a circular orbit of radius 0.086pc. The IMBH sinks toward the SMBH through dynamical 
friction, but dynamical friction becomes ineffective when the IMBH reached the radius inside which the 
initial stellar mass is comparable to the IMBH mass. This is because the IMBH kicks out the stars. 
This behavior is essentially the same as the loss-cone depletion observed in simulations of massive 
SMBH binaries. After the evolution through dynamical friction stalled, the eccentricity of the orbit 
of the IMBH goes up, resulting in the strong reduction in the merging timescale through gravitational 
wave radiation. Our result indicates that the IMBHs formed close to the galactic center can merge 
with the central SMBH in short time. The number of merging events detectable with DECIGO is 
estimated to be around 50 per year. Event rate for LISA would be similar or less, depending on the 
growth mode of IMBHs. 

Subject headings: black-holes: gravitational radiation: 



1. INTRODUCTION 

Recent obse r vation s l)Matsumoto et al\ I200H : 
iMatsushita et a,l\ 1200(1) suggested that intermediate- 
mass black holes (IMBHs) exist in some starburst 
galaxies. The first such object, M82 X-1 (t he br ightest 
source m Figure 1 of iMatsumoto et al\ (|2Q0!!)) lays 
200pc off the dynamical center of M82, and has the 
estimated minimum mass of 700 Mq. Several sce- 
narios have been proposed for the formation of this 
object iEbisuzaki et fl?]l2001t iM iller & Hamihon 2002; 
iKawakatu ll2002H Taniguchi et al. 2000). Both Ebisuzaki 
et al. and Miller and Hamilton argued that IMBHs are 
formed through stellar dynamical process and merging. 
The main difference between them is simply in the time 
at which the merging occurs. Ebisuzaki et al. assumed 
that most of merging process occurred while participants 
were main-sequence stars. On the other hand. Miller 
and Hamilton assumed that the IMBH grew through 
merging of smaller black holes. 

Which of the two processes actually occur depends 
mainly on the initial thermal relaxation time of the clus- 
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ter. iPortegies Zwart et al\ l^2004^ showed, using A^-body 
simulations, that the runaway merging of massive star 
occurs if the initial relaxation time of the star cluster 
is less than 4 Myrs, which is of the order of the life- 
tim e of massive sta.rs. T he compact star cluster MGG- 
11 (jMcCradv et a/.ll2003j) . which coincides with the lo- 
cation of M82 X-1, has the estimated relaxation time 
of 3 Myrs. On the other hand, the relaxation time of 
MGG9, which is more massive than MGGll, is signifi- 
cantly longer. This difference is consistent with the ex- 
istence and nonexistence of IMBHs in MGG-11 and 9, 
respectively. In the case of the scenario by Miller and 
Hamilton, the growth timescale would be much longer, 
and it is hard to explain why any IMBH can exist in a 
young cluster like MGG-11. 

Compared to the relaxation time of typical globu- 
lar clusters, which is 10^""^ yrs, the relaxation time 
of a few Myrs might sounds extremely short. How- 
ever, for young clusters, such short relaxation time 
is not unusual. For example, Arches and Quintuplet 
clusters, whose estimated ages are around 1-5 Myrs 
years, have the estimated re laxation time of 12M years 
' ijPortegies Zwart et aJ2002|) . Star clusters R136 in LMC 
and Westerlund 1 in our galaxy are other examples of 
such young compact clusters. 
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If the formation of IMBHs is not a rare event in 
young and compact clusters, two questions naturally 
arise. What is the final fate of the IMBH and its parent 
cluster, and whether or not it is related to the growth 
of the central black hole of the parent galaxy. The star 
cluster itself evolves through internal thermal relaxation, 
tidal truncation, and dynamical friction, much in the 
same way as globular clusters evolve. The main differ- 
ence is again that the timescale is much shorter. For 
example, if there is a cluster with mass IO^Mq at the 
distance of 30 pc from the center of our Galaxy, the dy- 
namical friction timescale would be 

2.38 X 10^ / r Y 
In A [m^cj 

l lOOkms-i H^lirj""' 

where u is a velocity dispersion and Mc is the mass of 
the cluster. 

As the cluster sinks toward the center of the parent 
galaxy, the tidal field of the galaxy becomes stronger. As 
a result, the cluster loses mass, and eventually becomes 
completely disrupted, leavin g the IMBH orbiting a round 
the central SMBH. IRS13E fMaiUard et a/."2004) looks 
like such a remnant cluster, comp osed of a single IMBH 
and several stars still bound to it ijPortegies Zwart et a,l\ 
12005ft IRS13E is located at 4 arcsec (~0.16pc) from Sgr 
A*. It appears as a cluster of seven individual stars 
within a projected diameter of '--^ 0.5 arcsec (0.02 pc). All 
these sources have a common westward proper motion, 
indicating that they are bound (at least six of them). The 
types of stars imply that it is a young star cluster with 
the age of a few Myrs. To keep these stars bound, the to- 
tal gravitational mass of IRS13E must exceed 13OOM0, 
about one order of magnitude larger than the estimated 
total mass of visible stars. One natural interpretation is 
that IRS13E is a remaining core of much more massive 
star cluster, with an IMBH of 1300Mo. 

In this paper, we consider the orbital evolution of 
an IMBH after its parent cluster is completely dis- 
rupted. The main question is the merging timescale 
of the IMBH and the central SMBH. In the case of 
a massive BH binary, which forms when two galax- 
ies each with a central massive BH merge, the merg- 
ing timescale has been th e area of active research since 
the p ioneering work by iBegleman. Blandford. fc ReesI 
l)1980D . In this cas e, the result of recent l arge-scale 
A^-body simulations ijMakino fc Funato1l2004j) suggests 
that the merging timescale is much longer than the 
Hubble time. They demonstrated that the evolu- 
tion timescale of the BH binary is proportional to 
the relax ation time of the parent ga l axy, as sug- 
gested by iBeglemanT Blandford. fc Reed l(T980|) . This 
conclusion is different from the results of previous 
simulations (Quinlan & Hernquist 1997; Makino 1993 
i Milosavlievic fc Merritt I I200lt l2003t IChatteriee et ah 
12003ft . but the difference is mostly due to the limi- 
tation in the number of par ticles in previous simular - 
tions. ISzell fc Merritt I 120051) and 'Berczik et al' ("2003) 
obtained similar results as Makino fc Funato ( 2004) , al- 
beit with somewhat smaller number of particles. 

However, we cannot directly apply these result to the 
case of an SMBH- IMBH binary, because of its very large 



mass ratio. In the case of SMBH-SMBH binary, the bi- 
nary need to interact with the field stars with the total 
mass comparable to that of the total mass of the binary 
to change its internal orbital parameters significantly. 
Thus, the loss-cone depletion occurs when the binary 
ejected out the mass comparable to its mass, and the 
structure outside the loss cone region is self-gravitating. 
In the case of SMBH- IMBH binary, the IMBH orbit can 
evolve by interacting with the mass comparable to the 
IMBH mass which is several orders of magnitude smaller 
than the SMBH mass. Thus, loss cone depletion can oc- 
cur when the IMBH ejected out the central stellar mass 
comparable to the IMBH mass, and the gravitational po- 
tential of this region, or actually of the region much fur- 
ther out, is dominated by the potential of the central BH. 
Thus, unlike the case of SMBH-SMBH binary, we need 
to study the evolution of a highly unequal mass binary, 
in the distribution of stars which are all bound to the 
primary component of the binary. 

The organization of the paper is as follows. In §2, we 
describe the numerical method and the initial conditions 
we used. In §3 and §4, we describe the result of the 
simulations. Summary and discussion are presented in 
§5. 

2. INITIAL MODELS AND NUMERICAL METHODS 
2.1. Initial Models 

The goal of this paper is to study the orbital evolution 
of an IMBH in the stellar distribution where the grav- 
itational potential of the central SMBH dominates. As 
the background stellar distribution, we adopt the stan- 
dard p cx r~''l^ cusp HBabcall fc Wolf M 9761 . One practi- 
cal problem with this distribution is that the total mass 
would be infinity if the density profile is given by this 
single power law. In addition, such a model is physically 
unacceptable, since the basic assumption for the Bahcall- 
Wolf cusp is that the gravitational potential is dominated 
by that of the central BH. 

In order to construct a model with central density slope 
of —7/4 and finite mass, we use Tre maine's t;- model with 
central BH (jTremaine et a/.l|l994j) . with one modifica- 
tion. The original Ty-model has the outer slope of —4. 
We constructed the model with outer slope of —5, just 
to make the mass in the outer region smaller. 

The density distribution of our modified model is given 

by 

ri M 

^''(-)-^ ,3-.(,g;,V2+i ^0<^?<3, (2) 

where is the total mass of field stars and tq is the 
scale length. The model with 77 = 5/4 correspond to 
profile p cx r and we use it in this paper. 

In order to construct an TV-body model in dynamical 
equilibrium, we need to construct the distribution func- 
tion (DF). We could obtain DF at least numerically by 
solving the Abel integral equation. However, since what 
we need is a model with an inner power-law cusp and 
some outer cutoff, it would be an overkill to obtain a dis- 
tribution function which exactly satisfies equation (|2Jl. 
So we approximate DF by the following formula. 

/W-/o(e)'/'(^g + e^)"'^ , (3) 
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where e 
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-l/(')+2) 



r;M^r(4 - 77) 



27/2^5/2^J-'Jr(5/2 - 77) 



and 



/i - 



8\/2 
7^ 



(4) 



(5) 



(6) 



DF of this form gives a correct asymptotic behavior for 
both e ^ 00 and e ^ 0. This formula has one adjustable 
parameter s. We chose s — 5 after some numerical tests. 

We generated the initial A^-body model in the follow- 
ing two steps. First, we generate positions of particles 
so that the density distribution obeys equation ^ . Then 
we assign the velocity to each particle, so that the veloc- 
ity distribution at the point of particle is consistent with 
the DF given in equation (|3J). Since the DF we used is 
an approximate solution, the constructed iV-body model 
is not in an exact dynamical equilibrium. However, as 
we will see in our initial model is in pretty good 

dynamical equilibrium. 

We chose the system of units in which the total mass 
of the system and gravitational constant are both unity 
and the total binding energy of the system is —1/4. We 
set the mass of SMBH to be Ms = 0.5, that of the IMBH 
Mj = 5 X 10"'*. The mass of the stellar distribution is 
therefore 0.4995. The stellar mass inside the radius r is 



(7) 



For the above choice of the system of units, rg ~ 2.39. 
When constructing the initial stellar distribution, we ex- 
clude the stars with periastron distance less than 10"'^, 
to avoid numerical problem. 

To convert the timescale obtained in our simulations to 
the physical timescale, we need to define conversions be- 
tween our system of units and real physical units. We use 
unit mass M = 6.0 x IO^Mq, unit length R ^ 0.86pc, 
resulting in the unit time of T = 4.6 x lO'^yrs. Mass 
and length units are chosen so that the stellar mass in- 
side the radius 0.7pc (assuming that the single power law 
with —7/4 slope con tinues to that radius) is 7.3 x IO^Mq 
(ICxenzel 61^ a/Tl200(l . Figure ^ shows the mass distri- 
bution for our mod el (assuming single power-law den- 
sity), the model bv iGenzel et al\ 12000), and observa- 
tional data. Note that our model is within the observa- 
tional error bars and prac tically indistinguishable with 
the model bv Genzel eTaP (2000). 

The IMBH particle is placed in a circular orbit at 
the radius 0.1 (0.086pc in physical units). Its mass is 
3000Mq. 

We performed four runs with different number of field 
particles (models Al to A4 of table 1). The mass of a 
field particle in model A4 is 30 M©, which is not too far 
from the mass of real stars (whatever they are) in this 
region. 

In table 1, models Bl to B4 were prepared to study 
the evolution of the IMBH after it approaches to less 
than O.Ol pc to the SMBH. For these models, we used rg 
much smaller than that is used for model Ax, to reduce 
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Fig. 1. — Cumulative mass distribution of Our Galaxy ob- 
tained from stellar and gas dynamics and that of our and Genzel 
2000 m odels. Triangles den ote mass estimates from the gas dy- 
namics IGuesten et aUll9871) . C rosses, circles and squares denote 
the isotropic mass modeling of |G«n^r et at\ fl99 6: 1997: 2000). 
Dashed curve is the mass model bv lGenzel et al\l200a) . whose stel- 
lar distribution is p(r) = po[l + (''Ao)^]~°/^, of po = 3.5x lO'^M©, 
ro = 0.17pc and a = 1.8. Solid curve is our model, drawn using 
single power-law mass distribution with p oc r~^/*. 



the total number of stars. We used ro = 0.1 and placed 
the IMBH at r = 0.01, so that the stellar distribution at 
the initial position of IMBH is still the power law with 
slope —7/4. The total stellar mass is 8.9 x 10^^, or about 
20 times the IMBH mass. Results of Bx models will be 
discussed in §4. Mass of a field particle in model B4 is 
3 Mq, which is of the same order with that of real stars 
in this region. Thus, two-body relaxation effect in this 
model is essentially the same as what would occur in the 
real galactic center. 

We made models CI an d C2 to test the effect of soft- 
ening parameters (see in ^2.2|) and Dl and D2 to test 
the validity of the initial model and relaxation effect (see 

2.2. Hardware and Numerical Integration Method 

For all calculations, we used simple direct-summation 
algorithm and fourth-order Hermit e scheme integrator 
with individual (block) time step ijMakino &: AarsethI 

For the calculation of gravitational forces from field 
particles (to both the field particle and black holes), 
we used the GRAPE-6 (and -6A) of Tokyo Univer- 
sity, a special purpo s e computer for iV — bodv simulation 
()Makino et adl2003l iFukushige et a/Jl2005j) . The calcu- 
lation of the forces from black holes was done on the host 
computer to maintain sufficient accuracy. We did not in- 
clude any relativistic effects and treated BH particle as 
massive Newtonian particle. 

In our calculation, we did not use regularization tech- 
nique to keep the calculation code simple. Thus, we need 
to apply softening parameters for gravitational interac- 
tions between particles. We use four different softening 
lengths, for BH-BH, SMBH-star, IMBH-star, and star- 
star interactions. 

For SMBH-IMBH interaction, we apply zero softening. 
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TABLE 1 

List of parameters for the simulations in this paper. 



Run 




Mi/m, 


n 




c 




{ 




g 


A 1 


yyyu 


inn 


n 1 


1 n— 
iU 


4 


1 n — 
iU 


6 


1 n— 


4 


A2 


19980 


20.0 


0.1 


10" 


4 


10" 


6 


10- 


4 


A3 


39960 


40.0 


0.1 


10" 


4 


10~ 


6 


10- 


4 


A A 
A4 


nnnnn 

yyyuu 


1 nn n 
iUU.U 


n 1 
U. 1 


1 n — 
lU 


4 


1 n — 
iU 


6 


1 n— 


4 


Bl 


1795 


100.0 


0.01 


10~ 


6 


10" 


8 


10" 


4 


B2 


3589 


200.0 


0.01 


10- 


6 


10- 


8 


10- 


4 


B3 


7177 


400.0 


0.01 


10" 


6 


10" 


8 


10" 


4 


B4 


17942 


1000.0 


0.01 


10" 


6 


10" 


8 


10" 


4 


CI 


99900 


100.0 


0.1 


10- 


4 


10- 


6 


10" 


3 


C2 


99900 


100.0 


0.1 


10- 


4 


10- 


6 


10- 


5 


Dl^ 


10000 






10" 


4 






10" 


4 


D2h 


100000 






10" 


4 






10- 


4 



system with a peak speed of 1 T-flops. For all calcula- 
tions, the total energy is conserved to better than 0.5% 
for all Ax, Cx and Dx runs. For models Bx the total 
energy conservation is shown in Figure [T51 

2.3. Stability and relaxation effect 
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They cannot easily come close enough to each other for 
the numerical difficulty to occur. If such close encounter 
occurs, the timescale of orbital evolution through grav- 
itational wave radiation becomes short enough so that 
our pure Newtonian treatment is not really valid. 

The softening for BH-star interaction should be deter- 
mined with some care, since close encounters do occur 
and too large softening can affect the orbital evolution of 
the IMBH. For the softening of SMBH-star interaction, 
we used es-s = 10"'' and 10~^ in Ax and Bx runs, re- 
spectively. These values are chosen so that the effect of 
softening is small enough for stars at the radius compa- 
rable to that of the IMBH. 

For IMBH-star interaction, we chose the softening so 
that it is smaller than 90-degree turnaround distance by 
at least two orders of magnitude at the initial condition. 
As the IMBH approaches to SMBH, the velocity disper- 
sion becomes larger, the difference between the softening 
length and the 90-degree turnaround distance becomes 
smaller, but it was always kept significantly larger than 
1. 

The criterion of the softening for star-star interaction 
is more complicated. Since the "stars" in our model is 
still significantly heavier than real stars (assuming we 
know what the mass of real stars in this region), it is 
desirable to use large softening to reduce the relaxation 
effect. On the other hand, the softening should be small 
enough not to affect the distribution of stars. For most 
of runs we used the softening length of 10^"'. As the test 
calculations, we performed two runs with e^-s = 10^"^ 
and 10^^, (runs CI and C2). As will be discussed in 
^2.'M these runs gave essentially the same results are the 
standard run (run A4). So we used €s-s = 10"'' for all 
other runs. Tabled also gives the softening parameters 
used. 

The largest calculation (model B4) took about five 
weeks on a single-host, single processor-board GRAPE-6 



Fig. 2. — Radial density profiles at times T = and 1000 for 
runs Dl and D2. 
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Fig. 3. — Lagrangian radii around SMBH for runs Dl and D2. 
From bottom to top, the radii enclose 5 X 10~^, 10~*, 2 X 10~*, 5 X 
10-*, 10-3,2 X 10-3,5 X 10-3, 0.01 and 0.02 of the unit mass. 



As described in §2.1 our initial model is not in exact 
dynamical equilibrium. In addition, the system would 
evolve through two-body relaxation even if the initial 
model is in exact dynamical equilibrium. To see these 
effects, we performed two test calculations (models Dl 
and D2), where we let the system without the IMBH to 
evolve for 1000 time units. Figures |21 and O show the 
result. 

We can see that the density profiles are practically un- 
changed, and that Lagrangian radii do not show any sys- 
tematic evolution. 

2.4. Effect of star-star softening 

As we discussed in §2.2, the choice of softening for star- 
star interaction might have some unpredictable effect on 
the evolution of the orbit of the IMBH. To see if there 
is any such effect, we performed two runs, models CI 
and C2, which started from the same initial condition as 
model A4 but with different values of es-s- 

The result is shown in figure 21 There is no systematic 
difference among these three runs. So we can conclude 
that the choice of tg-s has no significant effect on the 
result. 

3. RESULTS 

3.1. Hardening Rate 

Figure shows the time evolution of the semi-major 
axis of the IMBH, or the SMBH-IMBH binary Here and 
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Fig. 4. — Evolution of the semi-major axis and the eccentricities 
of the black hole binary. Dotted, dashed, and solid curves denote 
es-s = 10"'^, 10"^ and 10~*, respectively. 
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Fig. 5. — Evolution of the semi-major axis of the IMBH. Dotted, 
dash-dotted, dashed, and solid curves are the results of models Al - 
4. Thin curve shows the theoretical curve calculated equation 1161 . 



hereafter, we refer to orbital elements and other quan- 
tities of IMBH-SMBH binary as those of the IMBH, for 
simplicity. The semi-major axis of the IMBH is given by 

GMsMi 



(8) 



(9) 



2Eb ' 

where Ei, is the binding energy of the IMBH 

1 2 GMsMi 
E» = 2A^n 

Here, Vf, is the relative velocity of the two BHs and /i is 
the reduced mass defined as 

ti = MsMi/{Ms+Mi). (10) 

From figure |S1 we can see that the orbital evolution of 
the IMBH is practically independent of the number of 
stars in the parent galaxy. In all models, the IMBH is 
much more massive than the field stars. So this result is 
not surprising. 



The thin solid curve of in figure |5l is the theo- 
retical prediction for the evolution of the IMBH or- 
bit, obtaine d using standard dynamical friction formula 
ijBinnev fc T rcmainc 1987) 



dVb 
dt 



47r InKG^pMi 



G{X). 



(11) 



where 



and 



2X 



G{X) = erf (X) - —= exp(-X''), (12) 



X = Vf,/{V2a). 



(13) 



Here a is the velocity dispersion and we used the circu- 
lar velocity H = GMs /a as the velocity of the IMBH. 
We used InA = 8 here. For the calculation of dynami- 
cal friction in inhomogeneous background distribution, 
it has been suggested that taking the outer cutoff of 
Coulombs logarithm to the distance of the object from 
the ce nter of the parent stell ar system gives good esti- 
mate (iHashimoto et a?Jl2003j) . The lower cutoff is 90- 
degree turnaround distance, which is given by 



10^-^a. 



(14) 



Thus, we have log A ~ 7 independent of the location 
of the IMBH. Using equation 0, equation Ullfl can be 
rewritten as follows 



da 



10.71nAGi/VoM/a^/^ 



(15) 



where po is the stellar mass density at r = 1 and we 
assume the circular motion so that we used a = r. 

We chose X = 1 here. This differential equation has 
the analytic solution given by 



a = flo 



To-r 

To 



where 



To 



p. 1/4 ,^3/2 

0.37 Qq Mg 
In A GV2 p^Mi 



(16) 



(17) 



We can see that the agreement between the theoretical 
prediction and numerical result is pretty good, at least 
for the early period [T < 600). In the later phase, it 
seems the numerical results show the slowing down of 
the evolution. 

To see the difference between the theoretical predic- 
tion and numerical results more clearly, we calculated 
the hardening rate /S, defined as 



f3- 



Aa 
' At' 



(18) 



Here A(l/a) — \ai — ao|, where ai and ao are the semi- 
major axis of the IMBH at times t = to and to + At, 
respectively. We use At — 100 for all values of to- Fig- 
ure shows the result. The agreement between the the- 
oretical prediction and numerical results is fairly good 
for a > 0.01. For a < 0.01, numerical result gives the 
hardening rate smaller than the theoretical prediction, 
and the difference becomes larger as a becomes smaller. 

A natural explanation of this slowing down is the loss- 
cone depletion, similar to what happens in the case of 
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Fig. 6. — Hardening rate /3 in equation 1181 . plotted as a function 
of semi-major axis a at time (to + At/2), for r uns Al-4. Thin curve 
shows the theoretical prediction in equation 1151 . 
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Fig. 7. — Time evolution of Lagrangian radii of field stars. From 
bottom to top, the curves show the radii containing the mass 5 X 
10--^ (10-"' of the stellar mass), lO"*, 2 X 10""', 5 X 10"*, 10"^, 2 X 
IQ-"^, 5 X 10"^, 0.01 and 0.02. Four panels show the result of runs 
Al through A4 (model names are shown in panels). Thick solid 
curves show the semi-major axis of the IMBH. 



massive BH binaries. Figures |7| show the Lagrangian 
radii of field stars. We used the position of SMBH as the 
coordinate center. We can see that as the IMBH sink to- 
ward the center, the Lagrangian radius corresponding to 
the position of the IMBH starts to expand. For example, 
in the case of model A4, the radius enclosing the mass of 
2 X 10"'^ (fourth curve from the top) starts to expand at 
around t = 400, which is the time the IMBH semi-major 
axis crosses that radius. Radii enclosing smaller masses 
show similar tendency, though the expansion is faster for 
radii with smaller mass. 

Since the expansion of a Lagrangian radius starts when 
the IMBH reaches that radius, we can conclude that this 
expansion is due to the back reaction of dynamical fric- 
tion to the IMBH. Thus, when the stellar mass inside 



the IMBH semi-major axis becomes comparable to the 
IMBH mass, the effect of back reaction becomes signif- 
icant. Stellar mass inside r — 0.01 is about 0.1% of 
the total stellar mass, which is about the same as the 
IMBH mass. Thus, when the IMBH reaches the radius 
0.01, the effect of the IMBH to the stellar distribution 
becomes significant, and number density of field stars is 
reduced. This is the reason why the hardening rate be- 
comes smaller when a reaches 0.01. 

3.2. Change of the distribution of field stars 






Fig. 8. — Distribution of the field stars in the (J, E) plane at 
time T = for run A4 and those at T = 1000 for runs A4 and D2. 



Figures IHl show the distribution of the field stars in the 
(J, E) plane, for run A4 and D2 at T = and 1000. Here, 
E and J are the specific binding energy and specific total 
angular momentum of field stars. We take the center of 
mass of SMBH-IMBH binary as the origin tocalculate E 
and J. 

Note that we excluded field particles with the peri- 
astron distance to SMBH less than 10"'^ when we con- 
structed the initial condition. The left-hand-side cutoff 
in the distribution (at around J = 0.02) is due to this 
exclusion. 

When we compare the panels, it is clear that field par- 
ticles with small J ( J < 0.03) and large negative E {E < 
— 10) are depleted in model A4, while no such tendency is 
visible for model D2 (without IMBH). These particles are 
kicked out to high-energy, low-angular-momentum orbits 
by interaction with the IMBH. 

There are many particles in the area J < 0.1 and E > 
—0.2 in T = 1000 panel for run A4, while these are not 
in model D2. It is clear that these particles were kicked 
out by interaction with the IMBH from the small J and 
large negative E region. 

Figures El show the initial and final (T = 1000) total 
angular momentum of stars for runs A4 and D2. In the 
case of run A2, we can see that a fair number of particles 
which initially have small J (J < 0.2) got much larger 
J. It indicates that field particles are kicked out by the 
IMBH. The dispersion in the case of run D2 is purely 
due to the two-body relaxation. 
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SMBH binary eccentricity defined as 




10-' 



Fig. 9. — Distribution of the field stars in the (Jq, J(T)) planes for 
runs A4 and D2. Jq is the initial specific total angular momentum 
and J(T) is that at the time T = 1000. 




Fig. 10. — Distribution of the field stars in the {Eq, E(T)) planes 
for runs A4 and D2. Eq is the initial specific energy and E{T) is 
that at the time T = 1000. Top and bottom figure shows the 
distribution having a positive and a negative energy, respectively. 



Figures 1101 show the initial and final binding energies 
of particles for runs A4 and D2. Here, we can see the 
scatter is much larger for run A4. 

3.3. Evolution of eccentricity 




Fig. 11. — Evolution of the eccentricity of the IMBH for runs 
Ax. 



Figure shows the time evolution of the IMBH- 




L2 



{Ms + Mi)a 



(19) 



where L is the angular momentum of the binary. 

From these results, it is not clear if there is any system- 
atic change in e or whether or not the evolution depends 
on N . The fluctuation in e is bigger for small- TV calcula- 
tion (Al and A2 compared to A3 or A4). On the other 
hand, in run A4 e seems to show systematic increase after 
T = 600. As we can see from figure |S1 it may be related 
to the slowing down of the evolution of the semi-major 
axis, which is caused by the loss-cone depletion. 

In the next section, we will investigate this evolution 
of the IMBH orbit after the depletion of the loss cone. 

4. THE EVOLUTION OF IMBH ORBIT AFTER THE 
LOSS-CONE DEPLETION 

In the previous section, we have seen that the IMBH 
sink toward the SMBH through dynamical friction, and 
its orbital evolution slows down when the IMBH reaches 
the radius total stellar mass inside which is comparable 
to the IMBH mass. In other words, it occurs after the 
IMBH ejected out the neighboring stars. In this sec- 
tion, we analyze the evolution of the IMBH orbit after 
it kicked out the neighboring stars. In order to do so, 
we performed another set of simulations (run Bx), in 
which we placed the IMBH initially at the distance 1/10 
of that in runs Ax. Also, the stellar distribution is cut 
off at smaller radius, to reduce the total number of par- 
ticles. We performed runs Bl through B4. The mass of 
the field stars in run Bl is the same as that in A4. In 
B4, the mass ratio between the IMBH and field stars is 
10"^, and therefore the mass of field stars is around the 
solar mass. Two-body effects in run B4 is not larger than 
what would occur in real galaxies. 

4.1. Evolution of the semi-major axis 



10"' 





Bl 




B2 




B3 




B4 




theory 





















500 



1000 
T 



1500 



Fig. 12. — Evolution of the semi-major axis of the IMBH. Dotted, 
dash-dotted, dashed, and solid curves are the results of models Bl- 
4. T hin curve shows the theoretical curve calculated using equation 

Figure [T^ shows the evolution of semi-major axis. We 
can see that the evolution does not depend on the mass 
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of field stars. This result is not surprising since even for 
run Bl the initial relaxation time of field particles at the 
initial location of the IMBH is 7.3 x 10'' in the time unit 
and is much longer than the duration of the calculation. 
The slowing down is much more pronounced compared to 
that in runs Ax, simply because we followed the evolution 
of the IMBH down to smaller value of a. The evolution 
of the semi-major axis got stuck by T = 600, when the 
semi-major axis reached a ~ 0.003. 



0.01 




Figure [TH shows the evolution of the Lagrangian radii 
of field particles. As the IMBH sinks toward the SMBH, 
field stars expand. As the evolution of the IMBH slows 
down, the expansion of field stars also slows down. It is 
clear that the expansion, or depletion of the loss cone, is 
the cause of the slowing down of the IMBH orbit evolu- 
tion. 



10= 



-10^ -10^ -10° 



Fig. 15. — Distribution of the field stars in the (Eo,E{T)) plane 
for model B3 at T = 2000. Top and bottom panels show the 
particles with positive and negative energies, respectively. 



Fig. 13. — Energy error of the calculations. The curves give the 
results for runs Bx. 

Figure lOl shows the energy error for each of Bx runs. 
Errors are reasonably small before T — 500 for all runs. 
Quick increases in error are due to the increase in the ec- 
centricity of the IMBH. At T = 500, the binding energy 
of the IMBH is ~ 60% of the initial total binding energy 
of the system. Thus, relative energy error of 1% corre- 
sponds to the error of the semi-major axis of 1.6%, which 
is small compared to the overall change of the semi-major 
axis. 



10" 



10-' 



10-' 



10-- 



10-' 




500 1000 1500 2000 500 1000 1500 2000 

T T 

Fig. 14. — Lagrangian radii around SMBH. From bottom to 
top, the radii enclose 5 X IQ-^, 10"*, 2 X lO"' *, 5 X lO"-*, IQ-^, 2 X 
10-^,5 X 10"^, 0.01 and 0.02 of the unit mass. Thick solid lines 
indicate the semi-major axis of the IMBH. 



Figure El might give the impression that all field stars 
expand outward. Actually, that impression is wrong. 
Figure[T5lshows the initial and final energies of field stars 
for run B3. Large fraction of field stars show very small 
change in the energy, and some particles get very large 
energy. This is of course the same as what is visible in 
figure irni for run A4. Only the stars which strongly in- 
teracted with the IMBH are ejected, and other stars are 
essentially unaffected. The apparent expansion of outer 
Lagrangian radii in figure IT^ is the result of the ejection 
of the inner part of the field stars. 

4.2. Evolution of angular momentum 

Figure [TBI shows the evolution of eccentricities for runs 
Bl to B4. Unlike the case of runs Ax, here it seems 
clear that in all runs eccentricity shows slow increase in 
the early phase (T < 800), and in some runs there are 
phases when e approaches to very close to unity. How- 
ever, exactly when such phase occurs shows large run-to- 
run variation. 

This quick increase of the eccentricity is rather sur- 
pris ing, since in pr evious studies of m assive BH binaries 
^Makino & Funato 12 004: C hatteriee et ffli!.l2003D such in- 
crease has not been observed. However, since the config- 
uration of the system is quite different, evolution can be 
very different. In the case of SMBH binary, two massive 
binaries have similar mass, and field stars which interact 
with the binary are not strongly bound to the binary. 
However, in our case of SMBH-IMBH binary, all field 
stars are strongly bound to the SMBH, and have almost 
Keplarian orbits. Thus, celestial-mechanical effects such 
as the mean -motion resonance and Kozai mechanism 
llkozaillT96^ can play important roles. These effects, 
on average, work as the transport mechanism for angu- 
lar momentum from rapidly-rotating objects to slowly- 
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Fig. 16. — Evolution of the eccentricity of the black hole binary. 
The curves give the results for runs Bx. We plot (1 — e) for the 
vertical axis. 



rotating objects. In other words, these effects tend to 
increase the eccentricity of the IMBH. 




500 1000 1500 2000 500 1000 1500 2000 

T T 

Fig. 17. — Evolution of the normalized angular momentum of 
the BH binary. Dotted, dash-dotted, dashed, and solid curves give 
the results for L*,L*,L*, and L* , respectively. 

Figures 1171 show the normalized angular momentum 
of the IMBH, L* = L/Lq, where L is the angu- 
lar momentum calculated in Equation (jl9(l and Lq = 
y/ {Ms + Mj)a is the angular momentum of the circular 
orbit with the same semi-major axis. We can see that 
the slow increase of the eccentricity directly corresponds 
to the decrease of the total angular momentum, and the 
random walk of the angular momentum vector is small, 
especially for runs with large N. Apparently, only after 
the quick increase of the eccentricity the random walk 
of the angular momentum vector becomes noticeable. 
Thus, probably the high eccentricity state corresponds 
to some statistical equilibrium. 

4.3. Why the eccentricity goes up? 



Here we try to understand the mechanism which drives 
the increase of the eccentricity. We concentrate on run 
B4, since it has the largest number of field particles and 
therefore the statistical analysis of the behavior of field 
particles is the most reliable. In particular, as we can 
see in figure [T7I in run B4 the orbital plane of the IMBH 
remains close to the xy-plane, while in other runs the 
random walk of the angular momentum vector is signifi- 
cant. 

In order to understand the behavior of the IMBH, 
it would be useful to see with which field particles 
the IMBH interacted. Figure El shows the cumulative 
change of angular momentum of field particles as the 
function of their semi-major axis a. Here, the change in 
the angular momentum Al of a particle is defined as 

A/= [L(r)-L(T + Ar)].LiMBH , (20) 
where L(T) is its angular momentum vector at time T, 

and LiMBH is the unit vector with the direction of the an- 
gular momentum vector of the IMBH. We use AT = 50 
for all values of T. Thus this Al is the angular momen- 
tum change projected to the orbital plane of the IMBH. 
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Fig. 18.— Radial density profiles at times T = 0, 400, 800, and 
2000 for run B4. 

Figure El show the change of the density profiles for 
run B4. We can see that the central density decreases, 
and the inclination of the central cusp becomes flat. It 
is caused by the ejection of field stars by the IMBH. 

The difference between the plot for T = and T = 800 
is that, even though the semi-major axis of the IMBH is 
much smaller at T = 800, the field particles which ex- 
changed the angular momentum has much wider distri- 
bution at T = 800 than at T = 0. For T = 0, stars with 
the semi-major axis comparable to that of the IMBH 
dominate the change in the angular momentum, while 
for T — 800, stars with a comparable to that of the 
IMBH (0.003) have practically no contribution, and ef- 
fect of stars with a > 1 is not negligible. Of course, since 
the IMBH ejected most of stars with a < 0.01, it cannot 
interact with them. 

Figure 1201 shows the same cumulative plot but as the 
function of the periastron distance rp. For T = 0, par- 
ticles which go inside the IMBH orbit account for the 
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Fig. 19. — Cumulative change of angular momentum of field 
particles Al as a function of semi-major axis a for the times T = 
and 800. 




10"^ 10-' 10"- 10"' 10^10"^ 10"-' 10"- 10"' lo" 



Fig. 20. — Same as figure [T^ but as a function of the periastron 
distance rp for the times T = and 800. In the right hand panel, 
curves are, from top to bottom, total change, changes for particles 
with a <3.0, 1.0, and 0.3. 



change in the angular momentum. On the other hand, 
at T = 800 particles orbit outside the IMBH orbit are re- 
sponsible for the change in the angular momentum, and 
that tendency is more significant for field particles with 
small a. 

For field stars with relatively small a, we can draw the 
following picture. The IMBH has created the hole in the 
distribution function, and it can interact only with stars 
with semi-major axis larger or comparable to its apoc- 
enter distance. In other words, the IMBH can strongly 
interact with other stars only when it is at the apocen- 
ter. If it loses kinetic energy at its apocenter, it becomes 
strongly eccentric. 

Interaction with stars with large a is more complex. 
However, here the interaction can occur at positions 
other than the apocenter, since field particles with large 
a can have small rp. Thus, the increase of the eccentricity 
is likely to be driven by interactions with field particles 
with small a. In order to test this hypothesis, we per- 
formed one additional run, run E, which we started from 
the T = 800 snapshot of run B4 but we removed all stars 
with a > 0.3. Fig 1211 shows the result. The change in 
the eccentricity is largely similar. So we can conclude 
that the interactions with field stars with small a (and 
relatively large Vp) drive the increase of the eccentricity 
of the IMBH. 




T 



Fig. 21. — Comparison of run B4 and run E. As run E, we started 
from the T = 800 snapshot of run B4 but we removed all stars with 
a > 0.3. Left is the evolution of the eccentricity, and right is the 
evolution of the angular momentum of the BH binary. 



This mechanism of the incr ease of the eccentricit y 
is similar to that proposed bv iFukushige et al\ l)1992|) . 
They argued that the dynamical friction on an eccen- 
tric binary should be most effective at the apocenter 
and therefore an eccentric binary should become more 
eccentric. For the binary SMBH of comparable masses, 
this mechanism turned out to be ineffective, because the 
change in the orbital elements of a binary is really the 
result of rather complex three-body interaction between 
two binary components and the incoming star. However, 
in our case of IMBH-SMBH binary, the situation is very 
different. Since the mass of the IMBH is much smaller 
than that of SMBH, field stars must come close to the 
IMBH to change its orbit. Thus, the IMBH does have 
more chance to interact with field starts at apocenter 
than at pericenter, and its eccentricity increases. 

5. DISCUSSION AND SUMMARY 

5.1. Final fate of IMBH 

We have seen that the evolution of semi-major axis of 
the IMBH effectively stops, when it ejected field stars in 
that region. We have also seen that the eccentricity of the 
IMBH, after the evolution of semi-major axis stopped, 
can reach very close to unity. From the viewpoint of the 
evolution of SMBH, important question is what the final 
fate of the IMBH is. 




2 4 6 

[10*' yr] 



Fi g. 22. — Gravitational radiation timescale given by equation 
I23i for runs Bx. 
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Figure shows the timescale of the merging through 
gravitational wave radiation, for runs Bx. We converted 
the simulation units to physical units. 

The timescale of the merging through gravitational 
wave radiation is given by 



256 G^MsMi{Ms + Mi) 



F{e) (21) 



with F{e) is 



F[e)^{l~e^)' IH H 

^ ^ ^ ' I 24 96 



(22) 



ijPeters Ill964j) . The timescale tgr depends strongly on e, 
when e is close to unity. In our units, t^^ is given by 



V-6.3 X 10"F(e) 



O.Olpcy V3x106Mq 
1 

yr- 



Mr 



^sxiosmq; ^^^^ 

We can see from figure that the merging timescale 
can go down to less than 10^ yrs, and stays at that value 
for some time. In other words, if we take into account 
the effect of gravitational wave radiation, the IMBH will 
merge to SMBH. 

5.2. Event rate for gravitational wave detection 

iMatsubavashi et all l)2004f) estimated the event rate of 
the merging of IMBH-SMBH binary, assuming that all 
binaries eventually merge. Our present result indicates 
that this assumption of 100% merging is valid. They 
estimated the event rate to be 20 ~ 70 per year, for the 
detection limit of /i « 10~^^. Here, h is the dimensionless 
amplitude of gravitational wave. The frequency of the 
gravitational wave in their final merging phase is 1 0~^ to 
10^ H z. It is within the target range of LISA(LISA repOT^ 
I2OOOI) and DEC IGQiSeto et al. 2001). 

IMatsubavashi et aL (,200 4) considered two limiting 
cases for the growth of SMBH. In the first case, the 
growth is hierarchical. BHs always merge with another 
BH with a similar mass. In the second case, the growth 
is monopolistic and one BH always merge with IMBH 



of small mass. In the first case, large number of events 
has rather low amplitude, since they come from IMBH- 
IMBH mergings. In the second case, most events come 
from IMBH-SMBH merging, and therefore amplitude is 
bigger than that in the first case. DECIGO will be able 
to detect all events for both cases, while LISA would not 
detect the majority of events in the case of hierarchical 
growth. In the case of the monopolistic growth, both DE- 
CIGO and LISA wiU detect most of events. Thus, DE- 
CIGO event rate would be around 50, while LISA event 
rate would be 5-50, depending on the growth mode. 



5.3. Summary 

In this paper, we studied the orbital evolution of the 
IMBH after its parent cluster is completely disrupted. 
Our main findings are summarized as follows. 

Initially, the IMBH sink toward the SMBH through 
dynamical friction. However, the evolution of the semi- 
major axis stops when the IMBH approaches to the ra- 
dius at which the initial stellar mass is comparable to the 
IMBH mass. For our galaxy, it corresponds to around 
0.01 pc. If the IMBH remains in a circular orbit, the 
merging timescale by gravitational wave radiation would 
be 6.3 X 10^3 yrs. 

However, in this region the eccentricity of the IMBH 
approaches to unity, and therefore we expect the IMBH 
to quickly merge with the SMBH. This increase is due to 
interactions with field stars with periastron larger than 
the semi-major axis of the IMBH. The fact that many 
of the stars very close to the galactic center have large 
eccentricities is probably explained by the same mecha- 
nism. For example, S2 has e ^ 0.87 and a ~ 0.119"(~ 
4 X 10-^Dc;)( Schodc l et alSmi^ . The distribution of ec- 
centricities of known stars very close to SgrA* is not 
consistent with being isotropic. 
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